hide
# Set ggplot theme for figures
theme_set(theme_bw() +
theme(
# Titles
plot.title = element_text(hjust = 0, face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40"),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
# Axes
axis.line.y = element_line(),
axis.ticks.y = element_line(),
axis.line.x = element_line(),
axis.ticks.x = element_line(),
axis.text = element_text(size = 11),
# Facets
strip.text = element_text(color = "white",
face = "bold",
size = 11),
strip.background = element_rect(
color = "white",
fill = "#00736D",
size = 1,
linetype="solid"))
)
# Building a GMRI theme based on Wall street Journal and NYTimes theme
# base settings from {ggthemes}
theme_gmri <- function(base_size = 10,
bg_color = "lightblue",
base_family = "sans",
title_family = "sans",
facet_color = "teal") {
# Color from gmRi palette, sets background color
#colorhex <- gmRi::gmri_cols()[bg_color]
facet_hex <- gmri_cols()[facet_color]
# Set up theme
theme_foundation(
base_size = base_size,
base_family = base_family) +
theme(
# Major Elements
line = element_line(linetype = 1, colour = "black"),
rect = element_rect(fill = "transparent",
linetype = 0,
colour = NA),
text = element_text(colour = "black"),
title = element_text(family = title_family, size = 12),
# Axis elements
axis.text.x = element_text(colour = NULL),
axis.text.y = element_text(colour = NULL),
axis.ticks = element_line(colour = NULL),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = NULL),
axis.line = element_line(),
axis.line.y = element_blank(),
axis.text = element_text(size = 11),
# Legend Elements
legend.background = element_rect(),
legend.position = "top",
legend.direction = "horizontal",
legend.box = "vertical",
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
# Panel/Grid Setup
panel.grid = element_line(colour = NULL, linetype = 3, color = "gray70"),
panel.grid.major = element_line(colour = "black"),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
# Title and Caption Details
plot.title = element_text(hjust = 0, face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
plot.margin = unit(c(1, 1, 1, 1), "lines"),
# Facet Details
strip.text = element_text(color = "white", face = "bold", size = 11),
strip.background = element_rect(
color = "white",
fill = facet_hex,
size = 1,
linetype="solid"))
}
# Set the theme
# theme_set(theme_gmri(base_size = 10, base_family = "sans", title_family = "sans", facet_color = "teal"))
2021 Sea Surface Temperature
An Exceptional Year
2021 was an exceptionally warm year for the Gulf of Maine. It had the warmest fall on record for the area, and the second warmest summer. Much of the year the region experienced what would be considered “heatwave conditions”.
hide
# Load the timeseries
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)
# Clean up the data - add labels
region_timeseries <- region_timeseries %>%
distinct(time, .keep_all = T) %>%
mutate(time = as.Date(time),
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
week = lubridate::week(time),
doy = yday(time),
season = metR::season(month_num, lang = "en"),
#Set up correct year for winters, they carry across into next year
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
season_yr = ifelse( (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
)
# Get heatwave statuses for each day:
# Uses area weighted sst by default
region_hw <- pull_heatwave_events(temperature_timeseries = region_timeseries)
# Format dates and seasons and heatwave outcomes
region_hw <- region_hw %>%
mutate(time = as.Date(time),
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
week = lubridate::week(time),
doy = yday(time),
season = metR::season(month_num, lang = "en"),
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
#Set up correct year for winters, they carry across into next year
season_yr = ifelse(
(season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
)
# Set up heatwave status outcomes for fancy table features:
region_hw <- region_hw %>%
mutate(hw_outcomes = case_when(
mhw_event == TRUE ~ 1,
mcs_event == TRUE ~ 0,
TRUE ~ 0.5))
Departure from Historic Conditions
hide
# Pathto ERSST on box
ersst_path <- cs_path("okn", "ersst")
# Load temps and anoms (1982-2011 climatology)
ersst_temp <- stack(str_c(ersst_path, "ERSSTv5_merged.nc"))
ersst_anom <- stack(str_c(ersst_path, "ERSSTv5_anom.nc"))
# Crop GoM
# function to mask
mask_nc <- function(ras_obj, mask_shape){
# Get stats from ranks
m1 <- raster::mask(rotate(ras_obj), mask_shape)
m1 <- raster::crop(m1, mask_shape)
return(m1)}
# Load the shape, mask ERSST with it
gom_shape <- read_sf(poly_path)
gom_ersst <- mask_nc(ras_obj = ersst_temp, mask_shape = gom_shape)
gom_ersst_anom <- mask_nc(ras_obj = ersst_anom, mask_shape = gom_shape)
# Clean them up - temp
ersst_temp_tl <- as.data.frame(cellStats(gom_ersst, mean)) %>%
rownames_to_column(var = "Date") %>%
setNames(c("Date", "sst")) %>%
mutate(Date = str_sub(Date, 2, -1),
Date = str_replace_all(Date, "[.]", "-"),
Date = as.Date(Date),
yr = str_sub(Date, 1, 4))
# anomalies
ersst_anom_tl <- as.data.frame(cellStats(gom_ersst_anom, mean)) %>%
rownames_to_column(var = "Date") %>%
setNames(c("Date", "sst_anom")) %>%
mutate(Date = str_sub(Date, 2, -1),
Date = str_replace_all(Date, "[.]", "-"),
Date = as.Date(Date),
yr = str_sub(Date, 1, 4))
# Join them together
ersst_tl <- left_join(ersst_temp_tl, ersst_anom_tl) %>%
mutate(yr = as.numeric(yr))
# Get Yearly Average
ersst_yr <- ersst_tl %>%
group_by(yr) %>%
summarise(sst = mean(sst),
sst_anom = mean(sst_anom),
.groups = "drop") %>%
mutate(yr = as.numeric(yr)) %>%
filter(yr <= 2019)
# Same for OISST
oisst_yr <- region_hw %>%
filter(between(yr, 1982, 2021)) %>%
group_by(yr) %>%
summarise(sst = mean(sst), .groups = "drop")
# Make Splines
ersst_smooth <- as.data.frame(spline(ersst_yr$yr, ersst_yr$sst))
oisst_smooth <- as.data.frame(spline(oisst_yr$yr, oisst_yr$sst))
# Plot them
ggplot() +
geom_line(data = ersst_smooth, aes(x, y, color = "ERSSTv5"),
group = 1, alpha = 0.4, linetype = 1) +
geom_point(data = ersst_yr, aes(yr, sst, color = "ERSSTv5"),
size = 1) +
geom_line(data = oisst_smooth, aes(x, y, color = "OISSTv2"),
group = 1, alpha = 0.4, linetype = 1) +
geom_point(data = oisst_yr, aes(yr, sst, color = "OISSTv2"),
size = 1) +
scale_color_gmri() +
theme_gmri() +
labs(x = "", y = "Sea Surface Temperature",
color = "Temperature Data Record",
title = "Long-Term Temperature Record for the Gulf of Maine",
subtitle = "Current temperatures above 150-year highs") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
facet_zoom(xlim = c(1982, 2021)) +
theme(legend.position = "bottom")
Relationship to Climate Patterns
To better get a sense of what kind of year 2021 was it is helpful to relate it to what we would expect from similar years. One large factor in global weather patterns is the El Niño Southern Oscillation (ENSO).
hide
# Source: https://www.ncdc.noaa.gov/teleconnections/enso/soi
# library(rsoi)
# clim_idx <- download_soi()
soi_std <- tibble::tribble(#####
~"YEAR", ~"JAN", ~"FEB", ~"MAR", ~"APR", ~"MAY", ~"JUN", ~"JUL", ~"AUG", ~"SEP", ~"OCT", ~"NOV", ~"DEC",
1951, 1.5, 0.9, -0.1, -0.3, -0.7, 0.2, -1.0, -0.2, -1.1, -1.0, -0.8, -0.7,
1952, -0.9, -0.6, 0.5, -0.2, 0.8, 0.7, 0.5, 0.1, -0.2, 0.4, 0.0, -1.2,
1953, 0.3, -0.5, -0.2, 0.2, -1.7, 0.1, -0.0, -1.2, -1.2, 0.1, -0.3, -0.5,
1954, 0.7, -0.3, 0.3, 0.6, 0.5, 0.1, 0.4, 1.1, 0.2, 0.3, 0.1, 1.4,
1955, -0.5, 1.9, 0.6, -0.1, 1.0, 1.3, 1.6, 1.5, 1.3, 1.5, 1.2, 1.0,
1956, 1.3, 1.6, 1.3, 0.9, 1.4, 1.1, 1.1, 1.2, 0.1, 1.8, 0.2, 1.1,
1957, 0.6, -0.1, 0.2, 0.2, -0.7, 0.2, 0.2, -0.5, -0.9, 0.0, -1.0, -0.3,
1958, -1.9, -0.5, 0.3, 0.4, -0.5, 0.3, 0.4, 0.9, -0.3, 0.1, -0.4, -0.6,
1959, -0.9, -1.4, 1.3, 0.4, 0.5, -0.1, -0.3, -0.1, 0.0, 0.5, 0.9, 0.9,
1960, 0.1, 0.1, 1.0, 0.8, 0.5, 0.1, 0.5, 0.8, 0.7, 0.1, 0.5, 0.8,
1961, -0.3, 0.9, -1.8, 0.8, 0.3, 0.1, 0.2, 0.2, 0.1, -0.3, 0.5, 1.5,
1962, 2.0, -0.3, 0.1, 0.2, 1.1, 0.7, 0.1, 0.6, 0.4, 1.0, 0.3, 0.2,
1963, 1.0, 0.6, 1.1, 0.8, 0.4, -0.5, -0.1, 0.0, -0.6, -1.2, -0.8, -1.2,
1964, -0.4, 0.0, 1.1, 1.1, 0.2, 0.8, 0.6, 1.5, 1.3, 1.3, 0.2, -0.3,
1965, -0.4, 0.4, 0.8, -0.5, 0.2, -0.6, -1.8, -0.7, -1.3, -0.9, -1.5, 0.2,
1966, -1.3, -0.2, -0.9, -0.2, -0.4, 0.3, 0.1, 0.6, -0.2, -0.1, -0.0, -0.3,
1967, 1.7, 1.7, 1.2, -0.0, -0.0, 0.6, 0.2, 0.7, 0.5, 0.1, -0.4, -0.6,
1968, 0.5, 1.3, 0.1, 0.0, 1.2, 1.1, 0.7, 0.3, -0.3, -0.1, -0.3, 0.2,
1969, -1.5, -0.5, 0.4, -0.4, -0.2, 0.2, -0.5, -0.1, -1.0, -0.9, -0.1, 0.4,
1970, -1.1, -1.0, 0.6, -0.1, 0.4, 1.0, -0.4, 0.6, 1.2, 1.0, 1.6, 1.9,
1971, 0.4, 2.0, 2.3, 1.7, 0.9, 0.4, 0.2, 1.5, 1.4, 1.7, 0.5, 0.3,
1972, 0.5, 1.1, 0.6, -0.1, -1.6, -0.5, -1.4, -0.5, -1.4, -0.9, -0.3, -1.3,
1973, -0.3, -1.4, 0.7, 0.1, 0.4, 1.1, 0.6, 1.3, 1.2, 0.8, 2.6, 1.8,
1974, 2.4, 2.1, 2.4, 0.9, 1.0, 0.4, 1.1, 0.8, 1.1, 0.9, -0.1, 0.2,
1975, -0.5, 0.8, 1.6, 1.2, 0.6, 1.3, 1.9, 2.0, 2.1, 1.7, 1.2, 2.1,
1976, 1.4, 1.7, 1.7, 0.3, 0.4, 0.3, -0.9, -0.8, -1.1, 0.4, 0.7, -0.3,
1977, -0.4, 1.2, -0.5, -0.4, -0.5, -0.9, -1.1, -0.8, -0.8, -1.0, -1.3, -1.1,
1978, -0.3, -2.7, -0.2, -0.3, 1.4, 0.7, 0.6, 0.4, 0.1, -0.4, -0.0, -0.1,
1979, -0.4, 1.0, 0.1, -0.1, 0.5, 0.6, 1.3, -0.2, 0.1, -0.1, -0.4, -0.7,
1980, 0.4, 0.3, -0.4, -0.6, -0.0, -0.0, -0.0, 0.4, -0.5, 0.0, -0.3, -0.1,
1981, 0.4, -0.2, -1.3, -0.1, 0.8, 1.2, 0.8, 0.7, 0.3, -0.4, 0.2, 0.5,
1982, 1.2, 0.3, 0.6, 0.1, -0.3, -1.0, -1.5, -1.7, -1.7, -1.7, -2.6, -2.2,
1983, -3.5, -3.6, -2.4, -0.9, 0.6, 0.0, -0.6, 0.1, 0.9, 0.4, -0.1, 0.0,
1984, 0.2, 0.9, -0.2, 0.3, 0.2, -0.3, 0.2, 0.4, 0.1, -0.3, 0.3, -0.1,
1985, -0.3, 1.2, 0.8, 1.2, 0.4, -0.4, -0.1, 1.0, 0.0, -0.4, -0.2, 0.2,
1986, 1.0, -1.0, 0.5, 0.3, -0.2, 1.0, 0.3, -0.4, -0.5, 0.6, -1.2, -1.4,
1987, -0.7, -1.2, -1.3, -1.4, -1.3, -1.1, -1.4, -0.9, -1.0, -0.4, -0.0, -0.5,
1988, -0.1, -0.4, 0.6, 0.1, 0.9, 0.1, 1.0, 1.5, 1.8, 1.4, 1.7, 1.2,
1989, 1.5, 1.2, 1.1, 1.6, 1.2, 0.7, 0.9, -0.3, 0.5, 0.8, -0.2, -0.5,
1990, -0.1, -1.8, -0.4, 0.2, 1.2, 0.3, 0.5, -0.2, -0.7, 0.3, -0.5, -0.2,
1991, 0.6, 0.3, -0.7, -0.6, -1.0, -0.1, 0.0, -0.4, -1.5, -1.0, -0.7, -1.8,
1992, -2.9, -0.9, -2.0, -1.0, 0.3, -0.6, -0.6, 0.4, 0.1, -1.4, -0.7, -0.6,
1993, -0.9, -0.7, -0.5, -1.2, -0.3, -0.8, -0.8, -0.9, -0.7, -1.1, -0.1, 0.2,
1994, -0.1, 0.3, -0.7, -1.3, -0.7, -0.4, -1.3, -1.2, -1.6, -1.1, -0.6, -1.2,
1995, -0.4, -0.1, 0.8, -0.7, -0.4, 0.1, 0.4, 0.3, 0.3, 0.0, 0.0, -0.5,
1996, 1.0, 0.3, 1.1, 0.8, 0.3, 1.2, 0.7, 0.7, 0.6, 0.6, -0.1, 0.9,
1997, 0.5, 1.7, -0.4, -0.6, -1.3, -1.4, -0.8, -1.4, -1.4, -1.5, -1.2, -1.0,
1998, -2.7, -2.0, -2.4, -1.4, 0.3, 1.0, 1.2, 1.2, 1.0, 1.1, 1.0, 1.4,
1999, 1.8, 1.0, 1.3, 1.4, 0.2, 0.3, 0.5, 0.4, -0.1, 1.0, 1.0, 1.4,
2000, 0.7, 1.7, 1.3, 1.2, 0.4, -0.2, -0.2, 0.7, 0.9, 1.1, 1.8, 0.8,
2001, 1.0, 1.7, 0.9, 0.2, -0.5, 0.3, -0.2, -0.4, 0.2, -0.0, 0.7, -0.8,
2002, 0.4, 1.1, -0.2, -0.1, -0.8, -0.2, -0.5, -1.0, -0.6, -0.4, -0.5, -1.1,
2003, -0.2, -0.7, -0.3, -0.1, -0.3, -0.6, 0.3, 0.1, -0.1, 0.0, -0.3, 1.1,
2004, -1.3, 1.2, 0.4, -0.9, 1.0, -0.8, -0.5, -0.3, -0.3, -0.1, -0.7, -0.8,
2005, 0.3, -3.1, 0.3, -0.6, -0.8, 0.4, 0.2, -0.3, 0.4, 1.2, -0.2, -0.0,
2006, 1.7, 0.1, 1.8, 1.1, -0.5, -0.2, -0.6, -1.0, -0.6, -1.3, 0.1, -0.3,
2007, -0.8, -0.1, 0.2, -0.1, -0.1, 0.5, -0.3, 0.4, 0.2, 0.7, 0.9, 1.7,
2008, 1.8, 2.6, 1.4, 0.7, -0.1, 0.6, 0.3, 1.0, 1.2, 1.3, 1.3, 1.4,
2009, 1.1, 1.9, 0.4, 0.8, -0.1, 0.1, 0.2, -0.2, 0.3, -1.2, -0.6, -0.7,
2010, -1.1, -1.5, -0.7, 1.2, 0.9, 0.4, 1.8, 1.8, 2.2, 1.7, 1.3, 2.9,
2011, 2.3, 2.7, 2.5, 1.9, 0.4, 0.2, 1.0, 0.4, 1.0, 0.8, 1.1, 2.5,
2012, 1.1, 0.5, 0.7, -0.3, 0.0, -0.4, -0.0, -0.2, 0.2, 0.3, 0.3, -0.6,
2013, -0.1, -0.2, 1.5, 0.2, 0.8, 1.2, 0.8, 0.2, 0.3, -0.1, 0.7, 0.1,
2014, 1.4, 0.1, -0.9, 0.8, 0.5, 0.2, -0.2, -0.7, -0.7, -0.6, -0.9, -0.6,
2015, -0.8, 0.2, -0.7, -0.0, -0.7, -0.6, -1.1, -1.4, -1.6, -1.7, -0.5, -0.6,
2016, -2.2, -2.0, -0.1, -1.2, 0.4, 0.6, 0.4, 0.7, 1.2, -0.3, -0.1, 0.3,
2017, 0.2, -0.1, 0.9, -0.2, 0.3, -0.4, 0.8, 0.5, 0.6, 0.9, 0.9, -0.1,
2018, 1.1, -0.5, 1.5, 0.5, 0.4, -0.1, 0.2, -0.3, -0.9, 0.4, -0.1, 1.0,
2019, -0.0, -1.4, -0.3, 0.1, -0.4, -0.5, -0.4, -0.1, -1.2, -0.4, -0.8, -0.6,
2020, 0.2, -0.1, -0.1, 0.2, 0.4, -0.4, 0.4, 1.1, 0.9, 0.5, 0.7, 1.8,
2021, 1.9, 1.5, 0.4, 0.3, 0.5, 0.4, 1.4, 0.6, 0.8, 0.7, 1.0, 1.5
)
#####
# Pick a date centered in each month so you can plot them as one timeline
month_centers <- data.frame(Month = toupper(month.abb),
Month_num = str_pad(seq(1,12,1), width = 2, pad = "0", side = "left"))
# Reformat
soi_long <- soi_std %>%
pivot_longer(names_to = "Month", values_to = "soi_z", cols = JAN:DEC) %>%
left_join(month_centers, by = "Month") %>%
mutate(Date = as.Date(str_c(YEAR, "-", Month_num, "-15")),
ribbon_fill = ifelse(soi_z > 0,"above", "below"),
ribbon_above = ifelse(soi_z > 0, soi_z, NA),
above_ymin = ifelse(soi_z > 0, 0, NA),
ribbon_below = ifelse(soi_z < 0, soi_z, NA),
below_ymin = ifelse(soi_z < 0, 0, NA)
)
# The plotalllayers thing wasnt cutting it
# Plot the ribbon plot instead
cutpoints <- rev(seq(-4, 4, 1))
col_rdbu <- rev(RColorBrewer::brewer.pal(n = length(cutpoints), name = "RdBu"))
# fix GeomRibbon
GeomRibbon$handle_na <- function(data, params) { data }
# Make figure
ggplot(filter(soi_long, YEAR >= 1982), aes(x = Date, y = soi_z)) +
geom_ribbon(aes(ymin = above_ymin, ymax = ribbon_above, fill = "La Niña"), alpha = 0.8) +
geom_ribbon(aes(ymin = below_ymin, ymax = ribbon_below, fill = "El Niño"), alpha = 0.8) +
scale_fill_manual(values = c("La Niña" = col_rdbu[2], "El Niño" = col_rdbu[7])) +
geom_line(size = 0.25, color = "gray20") +
labs(x = "Date", y = "SOI", title = "ENSO Index",
caption = "Years overlapping OISST data displayed", fill = "") +
theme_gmri() +
theme(legend.position = "bottom")
When we look at average temperatures for the Gulf of Maine, and how that relates to the ENSO index, we found there was NO relationship. Meaning that a current month’s El niño/ La Nińa status did not correlate with warmer/hotter months.
A similar climate index to ENSO, with more relevance to our area is the North Atlantic Oscillation (NAO).
Looking at the NAO index, we found there was ___ relationship. Meaning that a current month’s El niño/ La Nińa status did ____ with warmer/hotter months.
Daily Temperature Record
hide
# Last Calendar Year
this_yr <- region_hw %>%
filter(year(time) == 2021)
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"MHW Threshold" = "coral3",
"Daily Climatology" = "gray30")
# Set the label with degree symbol
ylab <- "Sea Surface Temperature"
# Plot the last 365 days
hw_temp_p <- ggplot(this_yr, aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
# geom_textpath(aes(y = mhw_thresh, color = "MHW Threshold"), label = "Heatwave Threshold", hjust = 0.8, lty = 3) +
#geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
geom_textpath(aes(y = seas, color = "Daily Climatology"), label = "Climatology Mean", hjust = 0.5, lty = 2) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
theme(legend.title = element_blank(),
legend.position = "none") +
labs(x = "",
y = ylab)
# Plot the last 365 days - anomaly scale
linetype_key <- c(
"Sea Surface Temperature Anomaly" = 1,
"Heatwave Event" = 1,
"MHW Threshold" = 3,
"Daily Climatology" = 2)
# Same Plot as Anomalies
hw_anom_p <- this_yr %>%
mutate(
sst = sst,
seas = seas,
sst_anom = sst_anom,
mhw_thresh = mhw_thresh,
anom_thresh = mhw_thresh - seas,
anom_hwe = hwe - seas) %>%
ggplot(aes(x = time)) +
geom_segment(aes(x = time, xend = time,
y = 0, yend = sst_anom),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time,
y = anom_thresh, yend = anom_hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst_anom, color = "Sea Surface Temperature Anomaly")) +
geom_line(aes(y = anom_hwe, color = "Heatwave Event")) +
geom_line(aes(y = anom_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
#geom_line(aes(y = 0, color = "Daily Climatology"), lty = 2, size = .5) +
geom_textpath(aes(y = 0, color = "Daily Climatology"),
label = "Climatology Mean", hjust = 0.5, lty = 2, show.legend = FALSE) +
scale_color_manual(values = color_vals) +
scale_linetype_manual(values = linetype_key, guide = "none") +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
guides(color = guide_legend(override.aes = list(linetype = linetype_key), nrow = 2)) +
theme(legend.title = element_blank(),
legend.position = "top") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
labs(x = "",
y = "Temperature Anomaly",
caption = paste0("Climate reference period : 1982-2011"))
# Show figure
hw_temp_p / hw_anom_p
2021 Heatwave Events
hide
# Isolate 2021
hw_2021 <- filter(region_hw, yr == "2021")
# Get Event numbers
events_21 <- hw_2021 %>% drop_na(mhw_event_no) %>% distinct(mhw_event_no) %>% pull(mhw_event_no)
events_num <- length(events_21)
# Pull their data
events_dat <- filter(region_hw, mhw_event_no %in% as.character(events_21))
# How long on average
avg_duration <- round(nrow(events_dat) / events_num, digits = 0)
During 2021 there were three main marine heatwave events, including one that began on the new year. The average duration for each heatwave event was 75 days, or over three months long.
Heatwave Seasonality
hide
# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_hw %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# Polar plot comparing when each year experienced peak heatwave conditions
polar_dat <- grid_data %>% filter(year %in% c("2012", "2021")) %>%
mutate(hw_line = ifelse(mhw_event == TRUE, sst_anom, NA),
hw_mag = ifelse(year == "2012", hw_line * -1, hw_line))
# Ribbon?
polar_dat %>%
ggplot(aes(flat_date, y = hw_line, ymin = 0, ymax = hw_line)) +
geom_ribbon(fill = "darkred", alpha = 0.75) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
facet_wrap( ~ year) +
geomtextpath::coord_curvedpolar() +
# theme_gmri() +
labs(x = "", y = "Temperature Anomalies")
hide
# Side by side - need data separately
polar_dat %>%
ggplot(aes(x = hw_line, xend = 0,
y = flat_date, yend = flat_date, color = factor(year))) +
geom_segment(alpha = 0.4) +
scale_y_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0), ) +
facet_wrap( ~ year) +
scale_color_gmri() +
theme_gmri() +
# coord_flip() +
labs(x = "", y = "Temperature Anomalies", fill = "Year")
Following Each Heatwave
Do some breakdowns of how each heatwave looked: - where was heat concentrated - how long did it last - was anything unique about it?
Heatwave Frequency Heatmap
hide
# Prep the legend title
guide_lab <- "Sea Surface Temperature Anomaly \u00b0C"
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) * c(-1, 1)
# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
labs(x = "",
y = "",
"\nClimate reference period : 1982-2011") +
#scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent",
limit = limit) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme_classic() +
theme(legend.title = element_text(angle = 90))
# Assemble pieces
heatwave_heatmap
Temperature Horizons
hide
# make slight modification to grid_data
hori_data <- grid_data %>%
mutate(year = factor(year),
year = fct_rev(year))
# cut out outliers and set cutpoints
cutpoints <- hori_data %>%
mutate(
outlier = between(
sst_anom,
quantile(sst_anom, 0.25, na.rm=T)-
1.5*IQR(sst_anom, na.rm=T),
quantile(sst_anom, 0.75, na.rm=T)+
1.5*IQR(sst_anom, na.rm=T))) %>%
filter(outlier)
# Set origin
ori <- 0
# Manually pick cut points
sca <- seq(-4, 4, length.out = 9)
sca_bins <- str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "\u00b0C")
sca_labels <- rev(sca_bins)
# Function to plot one or more years
anom_horizon_plot <- function(hori_data, origin = ori, scale_cutpoints = sca, labels = sca_labels){
ggplot(hori_data) +
geom_horizon(aes(flat_date,
sst_anom,
fill = ..Cutpoints..),
origin = ori,
horizonscale = sca) +
scale_fill_hcl(palette = 'RdBu', reverse = T, labels = sca_labels) +
facet_grid(year~.) +
theme(
panel.grid = element_blank(),
panel.spacing.y=unit(0, "lines"),
strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
legend.position = "left"
) +
scale_x_date(expand=c(0,0),
date_breaks = "1 month",
date_labels = "%b") +
labs(x = '',
fill = "Temperature Anomaly \u00b0C")
}
# Plot Horizon Plot Using All Years
anom_horizon_plot(hori_data = filter(hori_data, yr <= 2021),
origin = ori,
scale_cutpoints = sca) +
labs(title = "Temperature Anomaly Horizons",
subtitle = "How deviation from the norm has become the new-normal",
caption = "Climatatological Reference Period: 1982-2011")
Comparing the Gulf’s Two Hottest Years
Early into 2021 it was apparent that the year was on-par with the previous title-holder for warmest year on record.
hide
# Reset the outliers
# cut out outliers and set cutpoints
cutpoints_2 <- hori_data %>%
filter(year %in% c("2012", "2021")) %>%
mutate(
outlier = between(
sst_anom,
quantile(sst_anom, 0.25, na.rm=T)-
1.5*IQR(sst_anom, na.rm=T),
quantile(sst_anom, 0.75, na.rm=T)+
1.5*IQR(sst_anom, na.rm=T))) %>%
filter(outlier)
# 2012
hori_2012 <- hori_data %>%
filter(year == "2012") %>%
anom_horizon_plot(scale_cutpoints = cutpoints_2) +
facet_wrap(~year, strip.position = "top")
# 2021
hori_2021 <- hori_data %>% filter(year == "2021") %>%
anom_horizon_plot(scale_cutpoints = cutpoints_2) +
facet_wrap(~year, strip.position = "top") +
theme(plot.caption = element_text(colour = "gray25", size = 6)) +
labs(caption = "Climatatological Reference Period: 1982-2011")
# Put them together as one figure
comparison_horizons <- hori_2012 / hori_2021 + plot_layout(guides = "collect")
# Do a bunch of formatting
comparison_horizons &
theme(
#Legend
legend.position = "top",
legend.key.height = unit(.5, "lines"),
legend.key.width = unit(3.75, "lines")) &
guides(fill = guide_legend(title = "Temperature Anomaly",
title.hjust = 0.5,
title.position = "top",
label.position = "bottom",
nrow = 1,
byrow = T, reverse = T)) &
plot_annotation(title = "Comparing the Gulf of Maine's Two Hottest Years",
subtitle = "Duration of temperature anomaly horizons, colored by their strength")
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))
Streamplot
hide
# Set Factor Order for Bins
bin_order <- c(
"Less than -1",
"-1 to -0.5",
"Within 0.5",
"0.5 to 1",
"Greater than 1")
# Set Colors
bin_colors <- RColorBrewer::brewer.pal(n = length(bin_order),
name = "RdBu")
bin_colors <- rev(bin_colors)
# Count Days in Bins
day_types <- region_hw %>%
mutate(
`Less than -1` = if_else(sst_anom <= -1, TRUE, FALSE),
`-1 to -0.5` = if_else(between(sst_anom, -1, -0.5), TRUE, FALSE),
`Within 0.5` = if_else(between(sst_anom, -0.5, 0.5), TRUE, FALSE),
`0.5 to 1` = if_else(between(sst_anom, 0.5, 1), TRUE, FALSE),
`Greater than 1` = if_else(sst_anom >= 1, TRUE, FALSE)) %>%
group_by(yr) %>%
summarise(
`Less than -1` = sum(`Less than -1`, na.rm = T),
`-1 to -0.5` = sum(`-1 to -0.5`, na.rm = T),
`Within 0.5` = sum(`Within 0.5`, na.rm = T),
`0.5 to 1` = sum(`0.5 to 1`, na.rm = T),
`Greater than 1` = sum(`Greater than 1`, na.rm = T),
.groups = "drop") %>%
pivot_longer(cols = "Less than -1":"Greater than 1", names_to = "Anomaly Strength", values_to = "day_total", values_drop_na = FALSE) %>%
mutate(day_total = if_else(is.na(day_total), 0L, day_total),
`Anomaly Strength` = factor(`Anomaly Strength`, levels = bin_order))
# Plot the Streamplot
hw_streamplot <- ggplot(day_types, aes(yr, y = day_total, fill = `Anomaly Strength`)) +
geom_stream(type = "proportion") +
#geom_segment(y = .5, yend = 0.5, linetype = 3, color = "gray50", x = 1981, xend = 2021, size = 0.25) +
scale_x_continuous(breaks = seq(1980,2030, by = 10),
minor_breaks = seq(1975, 2030, by = 5),
expand = c(.03,0)) +
scale_y_continuous(labels = scales::percent_format(), expand = c(0,0)) +
scale_fill_manual(values = setNames(bin_colors, bin_order)) +
labs(y = "", x = "", fill = "",
title = "A Warmer Gulf of Maine",
subtitle = "Fractions of Each Year Experienced as Hot/Cold Anomalies",
caption = "1981 not a complete years in the data*") +
theme_minimal(base_size = 10) +
theme(
axis.line.x = element_line(),
axis.ticks.x = element_line(),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
legend.position = "top",
legend.key.height = unit(.5, "lines"),
legend.key.width = unit(3.75, "lines"),
panel.grid.minor = element_blank(),
panel.grid = element_line(size = 0.2),
plot.margin = unit(c(.5,1,.3,.5), "cm"),
plot.title.position = "plot",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
legend.margin = margin(c(7,0,0,0)),
legend.justification = "center") +
guides(fill = guide_legend(title = "Anomaly Strength \u00b0C",
title.hjust = 0.5,
title.position = "left",
label.position = "bottom",
nrow = 1,
byrow = T))
# Show figure
hw_streamplot
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))
Characteristics of the Hottest Years
Multivariate figure comparing characteristics of both years: 2012 & 2021
hide
# Use all years so they can be ranked independently
# Metrics?
year_mets <- hori_data %>%
filter(year != "2022") %>%
group_by(year) %>%
summarise(
total_days = n(),
`Average Temperature` = round(mean(sst, na.rm = T), 2),
`Temperature Above Average` = round(mean(sst_anom, na.rm = T), 2),
`Cumulative Degrees Above Average` = round(sum(sst_anom, na.rm = T), 0),
hw_days = sum(mhw_event, na.rm = T),
hw_events = n_distinct(mhw_event_no),
peak_temp = max(sst),
peak_anom = max(sst_anom)) %>%
mutate(
hw_day_pct = round((hw_days/total_days) * 100, 2),
`Average Heatwave Length` = round(hw_days/hw_events, 1)
) %>%
pivot_longer(names_to = "Var", values_to = "Metric", cols = `Average Temperature`:`Average Heatwave Length`)
# Snipe out the top5 for each
year_met_ranks <- year_mets %>%
filter(Var %in% c("Average Temperature", "Temperature Above Average", "Average Heatwave Length", "Cumulative Degrees Above Average")) %>%
group_by(Var) %>%
slice_max(n = 5, order_by = Metric) %>%
ungroup() %>%
mutate(
yr_focus = ifelse(year == "2021", TRUE, FALSE),
year = reorder_within(year, Metric, Var))
year_met_ranks %>%
ggplot(aes(year, Metric, color = yr_focus)) +
geom_col(aes(fill = yr_focus), show.legend = F) +
geom_label(aes(label = Metric), show.legend = FALSE) +
facet_wrap(~Var, scales = "free", ncol = 2) +
coord_flip() +
scale_x_reordered() +
scale_y_continuous(expand = expansion(mult = c(0, .2))) +
labs(y = "",
x = NULL,
title = "How 2021 Stacks Up:",
subtitle = "Ranking Characteristics of Acute and Prolonged Warming Events") +
theme(plot.title = element_text(hjust = 0)) +
scale_color_gmri() +
scale_fill_gmri() +
theme(panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major. = element_line(linetype = 3, color = "gray50"))
Maps
- Where was 2021 the 1st, 2nd, third hottest year
Biological Events
- basking sharks
- hammerheads
- blue crab
A work by Adam A. Kemberling
Akemberling@gmri.org